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ABSTRACT 

We study the stability of rotating collisionless self-gravitating spherical systems 
by using high resolution N-body experiments on a Connection Machine CM-5. 

We added rotation to Ossipkov-Merritt (hereafter OM) anisotropic spher- 
ical systems by using two methods. The first method conserves the anisotropy 
of the distribution function defined in the OM algorithm. The second method 
distorts the systems in velocity space. We then show that the stability of sys- 
tems depends both on their anisotropy and on the value of the ratio between 
the total kinetic energy and the rotational kinetic energy. We also test the rele- 
vance of the stability parameters introduced by Perez et al. (1996) for the case 
of rotating systems. 
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1 INTRODUCTION 



Various analytical and numerical studies ([Antonov 1973| ; [Barnes et al. 19861 ; |Palmer fc Pa- 



[paloizou 1986| ; |Perez fc Aly 1996] ; Perez et al. 1996|) (and references therein) have shown that 
spherical, collisionless, self-gravitating anisotropic systems with components moving mainly 
on radial orbits are unstable. However, all these works considered non-rotating systems, and 
it is well known that rotation can play an important role in the dynamical evolution of sys- 
tems and modify their stability properties. It has been shown that rotation can be the cause 
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of the deformation of systems like globular clusters or weakly elliptical galaxies ( |Goodwin 
1997| ; [Staneva et al. 1996|) . 



The stability of rotating stellar systems is a very complex problem, Much work has 
been concerned with barred galaxies (which are rapidly rotating stellar systems) (for a 
review see, IAU Colloquium 157 ), but in a more general context, few studies have been 
devoted in the literature to this topic. For example, Papaloizou, Palmer and Allen (1991) 
have performed a series of numerical simulations to analyze the stability of systems where 
rotation was introduced by using the technique proposed by Lynden-Bell (1962). All their 
simulations produced endstates in which a triaxial bar appears. These important results 
cannot be considered general, since they were obtained for systems dominated by particles 
evolving on radial orbits, and was put in rotation by a specific procedure. In order to analyze 
the influence of rotation on the (in) stability of a given system, it is necessary to consider 
not only spherical systems with different kinds of anisotropy but also different methods for 
introducing the rotation. This paper develops such an analysis. We are also interested in 
testing the relevance of the stability parameters ( [Perez et al. 1996| ) on rotating systems. 



Perez et al. (1996) have shown that the stability of spherical self-gravitating non-rotating 
systems can be deduced from the 'anisotropic' component of the linear variation of the 
distribution function (see below Section 2.2). Such stability parameters can be computed 
from rotating systems. We show that they are still relevant for anisotropic systems as long 
as the rotational kinetic energy is not too large. 

The paper is arranged as follows. We describe in Section 2 the method that we use to 
obtain the initial non- rotating systems as well as the parameters describing the (in) stability 
of such systems. In Section 3, we detail the techniques used to introduce a parameterizable 
rotation to the initial conditions presented in Section 2. In Section 4, we show our numerical 
results on the (in) stability of various rotating systems generated with different procedures. 
Finally, the discussion and physical interpretation of our results are presented in Section 5. 

2 STABILITY AND INSTABILITY OF NON-ROTATING SYSTEMS 
2.1 Non-rotating Initial Conditions 



In a previous paper ( [Perez et al. 199£ ), we used the OM algorithm ( |Ussipkov 1979 ; Mcrritt 



|1985a| ; |Merritt 1985b| ; [Binney fc Tremaine 19871 ) to generate anisotropic self-gravitating 



spherical systems with various physical properties. This algorithm starts from a density 
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given by Pi SO ( r ) oc ip^ Q , where i(ji SO (r) is a known gravitational potential satisfying the Poisson 
equation, while n is the polytropic index (1/2 < n < 5). This density profile pi SO (r) is then 
deformed according to: 

Pani(r) := U + ^jjStoW, (!) 

where the anisotropic radius r a is a parameter which controls the deformation. 

Using the Abel inversion technique, this procedure allows one to define an anisotropic 
distribution function (hereafter DF) which depends both on E and L 2 through the variable 
Q := E + L 2 /2rl 

f (O) \/2 d r° alip iso dp ani 
JoW) 4tt 2 dQ Jq VV'iso - Q ' 

Once this DF has been computed, the initial conditions of our N-body numerical simula- 
tions are generated by choosing at random, from the above DF, the positions and velocities 
for the N particles. The density profile p an i( r ) defined by equation [I] is the probability den- 
sity from which the positions are generated. The velocities are generated from the velocity 
probability density deduced from the equation |2| (see appendix). 

It must be noted that, there is a fundamental limitation in the OM models: Any given 
value of the polytropic index n implies a critical value of r a below which the DF becomes neg- 
ative and unphysical in some region of phase space. Merritt (1985a) interprets this limitation 
as a simple illustration of the well-known fact that radial orbits cannot always reproduce an 
arbitrary spherical mass distribution. In theses cases, in order to extend the OM algorithm 
to highly radially anisotropic (r„ ~ 1) systems, we have arbitrarily set the DF equal to zero 
in this region. Such a procedure on DF affect only particles with a large value of Q. This 
procedure is applied for systems with a small value of r a which contains mainly particles with 
a small value of Q. Such a procedure affect then a very small number of particles (less than 
0.1% of the total number of particles). The systems with a modified DF are not strictly OM 
systems. However they conserve the properties which are for the present work: the density 
profils deduced from the modified DF are indistinguishable from the density profils given 
by equation (1) with the same value of r a , the Lindblad diagrams are very peaked around 
a small value of Q ([Perez et al. 199(j| ), they well correspond to highly radially anisotropic 
systems, and finally the radial dependence of the velocity anisotropy 
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is preserved. 

Finally, since each particle is initialized independently, the equilibrium DF f (E, L 2 ) of 
the system is in fact slightly perturbed. The perturbation is due to local Poissonian fluctua- 
tions of the density. The dynamical evolution of the system then represents the response of an 
anisotropic self-gravitating spherical equilibrium system submitted to such a perturbation. 

2.2 Stability Analysis 

The equilibrium DF of a spherical self-gravitating system depends only on the one-particle 
energy E and the squared total angular momentum L? . If gi denotes the perturbation 
generator, the linear variation of the DF can be written as 

5f = f§{^} + f§{^ 2 }- (4) 

If DF is a monotonic decreasing function^], the stability is then related to the Poisson 
brackets {gi, E} and {gi, L 2 } ( |Perez fc Aly 1996| ; |Perez et al. 1996|) . In our N-body simu- 
lations, these quantities appear as two random variables, e and A respectively, defined for 
each particle i ( [Perez et al. 199TS| ). The stability of the system can be predicted from the 
probability P t for e to be negative, and the statistical Pearson index ( palot 1973]) Pa of the 
variable A. All anisotropic collisionless self-gravitating non-rotating spherical systems with 

P A ^ 2.5 and P e Z 20% (5) 

are unstable, while those with 

P A ^ 2.5 and P e < 20% (6) 

are stable. The two other regions of the (P\P e ) plane, correspond to a transition between a 
stable system and an unstable system. In the particular case of OM models, these regions 
correspond to an anisotropy radius r a close to 2 (|Perez et al. 1996| ). 



3 GENERATION OF ROTATING SYSTEMS 



We consider all along this paper only systems with a DF which admits a monotonous decreasing dependence with respect 

Sfo ^ n and Ms. 



to all the isolating integrals of motion < 0, and ^& < 0) 
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3.1 Definition 

In order to generate virial-relaxed rotating spherical systems, we modify the non-rotating 
systems defined in the previous section by using techniques derived from the Lynden-Bell 
method (1962). Since this method preserves the position and the norm of the velocity of 
each particle, the systems are put in rotation without modifying their total potential and 
kinetic energy. In practice, we apply the following transformations to the velocity components 
{v r , Vg, v^} of the particles: 

Method 1 Method 2 
v r — > v r v r — > 0. 

ve — ► v e v e — > v e ( 7 ) 

V(f> ►l^l Vcf, - 



\ 



y2 



The first method then conserves the radial anisotropy defined in the OM algorithm, while 
the second method distorts the system in velocity space. 

The amount of rotation introduced by these methods can be evaluated through the ratio: 

fj, = K rot /K tot , (8) 

where K to t is the total kinetic energy, and K rot is the rotation kinetic energy defined by 
Navarro and White (1993): 

2 i= i [rf - (r< • L tot ) 2 ] 
Here, Lj is the specific angular momentum of particle i, L to t is a unit vector in the direction 
of the total angular momentum of the system, when the system does not rotate at all the 
Ltot vector is the null vector, and N is the total number of particles. In order to exclude 
counter-rotating particles, the sum in equation ([|) is actually carried out only over those 
particles satisfying the condition (Lj • L to t) > 0. 

In order to have systems with different strengths of homogeneous rotation (HR), we have 
applied either Method 1 or Method 2 to a fraction r of the total number of particles. This 
fraction has been constructed by choosing the particles at random in the overall system. 
When t — > 0, the system does not turn while, when r — > 1, the system rotation reaches its 
maximum value. 

In principle, there is no reason to consider only the case of homogeneous rotation. More- 
over, in order to roughly model the presence of a rotating massive object like those sometimes 
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Figure 1. n vs r for n = 4; HRi,r a = 1.5 (thin solid line); HRi, r a = 100 (thin dashed line) and HR2,r a = 1.5 (bold Solid 
line); HR 2 , r a = 100 (bold dashed line) 

considered in the center of some elliptical galaxies, we also consider inhomogeneous rotations 
(IR). In this case, we apply the above velocity transformations only to those particles placed 
at a radial distance smaller than k x Ri , where A; is a positive parameter and R 1 is the radius 
containing half of the system mass. If k — 0, the system does not turn while, if k — > +00 
the rotation has its maximum value. We have then four possible procedures to introduce 
a rotation motion on our initial conditions. The first two possibilities introduce a homoge- 
neous rotation by choosing particles at random in the whole system and modifying their 
velocities according either to the method 1 or 2, what defines the HRi and HR2 procedures, 
respectively. The other two possibilities introduce an inhomogeneous rotation by applying 
either the method 1 or 2 to modify the velocities of those particles placed within a given 
radial distance, which defines the procedures IRi and IR 2 , respectively. 

Figures [1] and ^ show the values of \i (equation ||) obtained from the four possible 
procedures. As we can see from such figures, only the HR 2 and IR 2 procedures lead to large 
fractions (/i > 10%) of kinetic rotation energy. We also note from these figures that the 
dependence of /x on r a depends on whether velocities have been modified according to the 
method 1 or 2. As a matter of fact, for the HRi and IRi procedures, the amount of rotation 
obtained (/i) is greater for large r a values than for small ones. On the contrary, for HR 2 and 
IR 2 , \i is larger for small r a values than for large ones. 
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Figure 2. fi vs k for n = 4, IRi, r a = 1.5 (thin solid line); IRi, r a = 100 (thin dashed line) and IR2, r & = 1-5 (bold Solid 
line); IR 2 , r a = 100 (bold dashed line) 

4 INFLUENCE OF THE ROTATION ON THE STABILITY 

Using the N-body code described in Alimi and Scholl (1993), we have performed on Connection- 
Machine 5 a series of numerical simulations^ of the evolution of the systems defined in the 
previous section. As the collisionless hypothesis is fundamental for interpreting our results, 
we have not continued our simulations beyond a few hundred dynamical times in order to 
avoid the later evolution where two-body relaxation arises. However, all our models reach 
a steady state before about 50 (where the initial dynamical time is estimated by the 
following formula = \/Yj r i I 'J2 v i the summations on initial positions and velocities are 



done on all the particles). We will then present our results for this interval. 

The physical mechanism of the radial-orbit instability for collisionless self-gravitating 
systems is well known. It has been described in detail by several authors (see (palmer 1994|) 
for example). The morphological deformation of the initial gravitational system resulting 
from this instability is mainly due to the trapping of particles with a low angular momentum 
in a bounded area of space. This trapping favors a deformation of the initial spherical 
system to an ellipsoidal or even a bar-like structure. To evaluate such a deformation, it is 



convenient to use the axial ratio defined from the moment of inertia tensor / (|Allen et al. 



1990|) . From the three real eigenvalues of J, Ai > A2 > A3, we compute the axial ratios 
O'l = A 2 /Ai and a 2 = A 2 /A 3 . These two quantities, which can always be defined because 
these eigenvalues never vanish, satisfy a\ < 1 < a 2 . In order to discriminate clearly between 
a bar-like structure, a quasi-sphere and a disk-like structure , we define the quantity / from 
a\ and a 2 

t The set of numerical simulations performed have been made with 16384 particles. Some experiments have been performed 
using more particles (65536), no significant change in comparison with the work presented here have been obtained. 
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Table 1. Stability parameter evolution for the first method for rotating systems with n = 4 



A bar-like structure is characterized by a,\ < 1 and 02 — 1 which implies a / value signifi- 
cantly larger than 1. A disc-like structure is characterized by a% ~ 1, a-i > 1 and < / < 1. 
Any system with a / value of order unity has a quasi-spherical structure. 



4.1 Rotating systems according to Method 1 

This type of rotation preserves the anisotropy of the non-rotating OM systems. The dis- 
tribution function of the rotating system depends only on the variable Q = E + L 2 /2r% 
as in the case of the non-rotating systems. We see in Table 1 that, in the case n = 4, the 
stability parameters defined in Section 2.2 are very weakly modified whatever the r and k 
parameters values are, that is to say, whatever the rotational kinetic energy is (low with this 
method). According to the conditions given by equations ||] and || we expect the (in)stability 
of the original non-rotating systems not to be modified when they are put in rotation. Our 
numerical simulations confirm this. In figure |3], we see that the evolution of axial ratios is 
similar for the rotating (dashed and dotted lines) and non-rotating (solid lines) systems. 
This results holds for the homogeneous and inhomogeneous rotations and whatever the n 
parameter value is. 

4.2 Rotating systems according to the second method 

The situation is now more complicated because the rotation procedure modifies the system's 
anisotropy. In the second method, the stationary OM distribution function is modified by 
a positive definite and time-independant transformation. The resulting DF is then always 
stationary and positive definite. Moreover as the modified systems are always spherical (no 
modification on positions have been performed), the new DF depends only on isolating 
integrals of motion, the energy and the squared angular momentum QPerez fc Aly 19961) . 



Figure 4. The evolution of the Virial ratio for the initial systems defined respectively by n = 4, r a = 1.5, and n = 4, r a = 100 
which have acquired their rotation according to the procedure HR2 with r = 30% (thin line) and IR2 with k = 1.4 (bold line) 
. A similar evolution for the Virial ratio (remaining very close to 1) is obtained for all runs. 



It is therefore a stationary solution of the collisionless Boltzmann-Poisson system. We have 
also verified that the resulting systems are always virialized, as confirmed by figure |4]. 

Let us first consider rotating systems with high values of r a . We recall that non-rotating 
OM systems with the same anisotropy parameter are stable ( Perez et al. 1996|) . The dy- 
namical evolution obtained for such systems in our present simulations (see Figure [|, top 
panels) allows us to distinguish two classes of behavior. When r a is large and the rate of 
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rotation stays modest (typically fi < 10%), we find that systems remain stable and spherical 
(ai ~ a2, / — 1). However, when r a is large and the rate of rotation becomes important 
(typically /i > 10%), we find that initially spherical systems develop a very soft bar-like 
instability (a 2 — 1, ai — 0.85, / ~ 1.3). 

Systems with small r a values have instead a different behavior, which depends on whether 
the rotational motion has been introduced by using a homogeneous or an inhomogeneous 
procedure. In the first case (HR 2 procedure), Figure | (bottom-left panel) shows that systems 
which are radial orbit-unstable without rotation (e.g., a\ ~ 0.55, a 2 ~ 1, / ^> 1) (solid line), 
become quasi-spheroidal (ai ~ 0.85, a 2 — 1.1, / — 1.3) when they have a modest rotation 
motion (/x < 10%) (dotted line). However, when rotation is important (typically /i > 10%), 
such systems develop a disc- like instability (ai — 1, a 2 — 1.25, / — 0) (dashed line). In the 
second case (IR 2 procedure) (bottom-right panel), the fact that rotation has been introduced 
only in a central region of the system prevents one from obtaining quasi-spherical systems 
and, therefore, the radial-orbit instability persists (ai — 0.65, a 2 ~ 1.1, / ~ 3.5) for systems 
with a modest amount of rotation (n = 4;IR 2 ;r a = 1.5;fc = 0.6) (dotted line). When the 
amount of rotation is high enough (// > 10%) a disc-like structure appears (a x ~ 1, a 2 ~ 1.25, 
/ ~ 0), as in the HR 2 procedure. The evolution of the axial ratio for the system (n = 4; 
IR 2 ; r a = 1.5; k = 1.4) is represented by dashed line. These results hold whatever the n 
parameter value is. In practice we have performed numerical simulations for three values of 
n(n = 3.5,4,4.5). 

Are the stability parameters P e and Pa still discriminating for such rotating systems?. We 
can see in Table 2 that non-rotating stable systems (r a = 100) are predicted to remain stable 
when the second rotation method is applied (P e stays smaller than 20 and Pa stays larger 
than 2.5) whatever the r and k parameters values are. On the other hand, non-rotating 
radial-orbit unstable systems (r a = 1.5) are predicted to become stable when sufficient 
Method 2 rotation is applied, from r = 0.4 and k = 0.8 which correspond on figures [I] 
and |2| to typically fi ~ 5%. As a matter of fact, in this case, P e becomes less than 20 
and P\ becomes greater than 2.5. Consequently, these stability parameters which have been 
constructed to predict the stability of non-rotating spherical systems are not relevant as 
soon as the quantity of kinetic rotation energy becomes large, typically \i ^ 10%. 
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Figure 5. The axial ratios vs. dynamical time for n = 4; r a = IOO-HR2 (top-left panel), r a = IOO-IR2 (top-right panel), 
r a = I.5-HR2 (bottom-left panel) and r a = I.5-IR2 (bottom-right panel). The amount of rotation is represented by using 
different kinds of lines: solid lines : r = 0(HR) or k = (IR), dotted lines: r = 0.3(HR) or k = 0.6 (IR), and dashed lines: 
t = 0.8(HR) or k = 1.4 (IR). 
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Table 2. Evolution of the stability parameter for systems with n = 4 and put in rotation by using the second method 

5 PHYSICAL INTERPRETATION AND CONCLUSIONS 

The rotational properties of collapsed systems depend to a large extent on the amount of 
angular moment before the collapse. In order to study in a realistic way the importance of 
rotation for the dynamics of self-gravitating systems, it is necessary either to attempt an 
analytical approach, or to perform a complete numerical study modelling the collapse and 
relaxation phases prior to the two-body relaxation phase. However, although the collapse 
of a system can be studied by using the introduced amount of rotational kinetic energy 
as parameter, it is difficult to extract general conclusions from this kind of experiments. 
As a matter of fact in this way the post-collapse physical features of the object cannot be 
completely controlled and hence it can then be difficult to study with these methods the 
influence of the rotation on post-collapse systems. 
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This justifies the method that we have used in this paper. Starting from virialized sys- 
tems with exactly known dynamical properties, we can study the influence of rotation by 
controlling its features. If the initial systems cover a wide variety of physical properties, 
and if the methods to introduce the rotation preserve certain fundamental features of these 
systems (invariance of mean energy, conservation or controlled modification of the distri- 
bution function), the numerical study will then be able to be used as a model to extract 
some general conclusions. As a matter of fact, our simulations start from a wide variety of 
initial conditions fully controllable through the dependence on the two parameters n and r a . 
On the other hand, the techniques used to introduce rotation to the systems preserve (as 
explained in Section 3) the properties that ensure that our models are spherical stationary 
solutions of the collisionless Boltzmann-Poisson system. 

The main properties found in our study are the following : 

• There do not exist spherical self-gravitating systems in "fast" rotation. Our simulations 
show in fact that, when fi 10%, the systems do not remain spherical but become lengthened 
along one or two axes depending on whether they are isotropic or anisotropic, respectively, 
when they do not have a rotational motion. 

• Rotation (in our case HR 2 and IR2) can allow for a reorganization of systems in velocity 
space able to modify their dynamical behavior. We have in fact shown that a moderate 
rotation (typically < /i ^ 10%) can stabilize and confer a quasi-spherical structure to 
systems that, when they are not rotating, suffer a radial-orbit instability. Therefore, there 
can exist rotating spherical self-gravitating systems. This is the case for our models with 
r a 2 and ji ^ 10%. 

• We have finally found that the stability parameters introduced in ( Perez et al. 1996 ) 
remain valid as long as /i ^ 10%. 

6 APPENDIX: GENERATION OF INITIAL CONDITIONS 
6.1 The initial positions for the particles 

Let us consider a density p iso given by the polytropic relation p iso = c n ip^ so , where 





and ifji SO is the solution of the Lame-Emden equation 
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This isotropic model is then deformed according to 

Pani{r) = (1 + —)p iso 
'"a 

where the anisotropic radius r a is a parameter which controls the deformation. The polytropic 
index n is chosen in the range ]0.5, 5] in order to the system admits a finite mass M (< r) = 
47r Jq r 2 p an i{r)dr . The total mass of the system is normalized to unity and we then compute 
for a large set of particles (1 < i < N) from the inverse function of M(x) and sin(x), the 
components 

U = r max M~ l (x) 
6i = 2 arcsin(\/x) 
0i = 2n.x 

where x is an uniform random variable on [0, 1]. The size of the system r max is chosen such 
that p ani {r > r max ) < 10~ 5 . 

6.2 The initial velocities for particles 

Let us consider the velocity components in spherical coordinates (v r , vg, v$), we compute 



ft = yv 2 + v 2 and a = arctan(fe/f ^). The probability p(r)for finding a particle in a 
volume <ir := dr d9 d<p dv r dvg da of the phase space is defined from the DF of the system 

p(T)dT = —f(T)r 2 dr sin 9d9d(f>v t dv t dv r da (11) 

In the Ossipkov-Merritt model DF is a function only on Q variable 

f = f(Q) with Q = \vl+ l (l + r l\ v 2 + ^r) (12) 



the equation (|TT| ) then reduces to 

piT)dT = J 8?I " ^ — f(r, v r , v t )dr dv t dv r (13) 
r, v r and v t are dependant random variables. In order to continue the integration of equation 



(|T3|) we introduce the variables R and (3 defined as following 

v r = R. cos (3 



^2 

1 + — = R. sin /3 

^' a 



We then get 



8nr 2 v t IR 2 , . . \ SvrVi? 2 sin/3 



f(r,v r ,v t ) — = f I — + V(r) I -p ^ydrdRdp 

N \ 1 + 7i) 
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We see from the previous expression that the random variable (3 is r and R independant 
and we have 

p(P)dp = ^fdp 

ic 2 2n2 / d2 \ 

p(r, i2)drd/2 = y-/ — + ^(r) drdi? 

' a 

The conditional probablity for finding a particle with a velocity defined by R at a given 
distance r is then 



| ro) = ^ ,°' - with p(r ) ' 



p(r ) 



and finally 

P(R \r = r ) 



Ait 
2tt 



M(= 1) 



(14) 



7— prr/ R' 2 f(^ + ^(r ))dR' 



(15) 



p(r G ) (l + |) ^(r-) y/2(Q-ijj(r )) 
We are now able to assign a velocity for each particle which the position have been previously 
determined. 

Ri = P (x\r = r a ) 
Pi = 2sin- 1 ( v /^) 
aci = 2tt.x 

where x is an uniform random variable on [0,1], and P~ l is the inverse function of the 
probability defined by equation (\TB±). Finally 

V ri = Ri.COSf3i 
V8i 



v<t>., 



\ 



\ 



1 + -f.i2j. sin/3j cosctj 



1 H — sin Pi sin aij 
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